!> this program solves K*phi = S where S is nuSigF*phi/k or Sext
!
! ./diff_2d biblisinp -ksp_type gmres -ksp_gmres_restart 10
!                      -pc_factor_levels 3 -pc_factor_reuse_ordering -pc_factor_fill 3
program diff_2d

    use class_input, only: indat, readinput
    use class_dat, only: rhsdat,coefdat,flux
    use class_2d_matx
    use class_solver
    use class_output

    implicit none
#include "finclude/petsc.h90"

    !> initiate input
    type(indat)  :: input
    !> initiate output (phi has data type and methods)
    !type(outdat) :: phi
    !> array of matrices for multigroup diffusion operator
    Mat, dimension(:), pointer   :: K
    !> initiate rhs stuff (includes fission, chi, and absorption)
    type(rhsdat) :: rhs
    !> absorption and other auxilliary data needed for rf's
    type(coefdat) :: dat
    !> flux
    type(flux)    :: phi
    PetscInt      :: ierr
    PetscMPIInt   :: iam, num
    double precision :: t1,t2
    call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

    t1 = MPI_Wtime()

    call MPI_Comm_size(PETSC_COMM_WORLD,num,ierr)
    if (num .ne. 1) then
       call MPI_Comm_rank(PETSC_COMM_WORLD,iam,ierr)
       if (iam .eq. 0) then
          write(6,*) 'This is a uniprocessor code only!  Goodbye.'
       endif
       SETERRQ(1,' ',ierr)
    endif

    !> read the input
    call readinput(input)
    !> build the problem matrix
    call buildk(K,input,rhs,dat)

    call solver(input,K,phi,rhs)
    call output(input,phi)
    ! call dealloc
    !write(0,*) "phi = "
    !write(0,'(4f10.5)') phi%val
    t2 = MPI_Wtime()

    write (*,*) " *** total time = ",t2-t1," seconds *** "

    call PetscFinalize(ierr)

end program diff_2d
